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In recent years it has been found that quantum systems can posses entanglement in equilibrium 
thermal states provided temperature is low enough. In the present work we explore a possibility 
of having entanglement in nonequilibrium stationary states. We show analytically that in a simple 
one-dimensional spin chain there is entanglement even at the highest attainable energies, that is, 
starting from an equilibrium state at infinite temperature, a sufficiently strong driving can induce 
entanglement, even in the thermodynamic limit. We also show that dissipative dephasing, on the 
other hand, destroys entanglement. 



I. INTRODUCTION 

It has been realized that entanglement is rather ubiq- 
uitous at sufficiently low temperatures. Ground states 
of interacting quantum systems typically posses entan- 
glement between different parts of a system. Similarly, 
when a system is at thermal equilibrium at sufficiently 
low temperature entanglement is often present. This has 
been studied in a number of works, see e.g., Refs. [1] and 
references in [2], and is called thermal entanglement [3]. 

Much less in known on the other hand about entan- 
glement in nonequilibrium stationary states. The main 
difficulty is that solving a nonequilibrium system, cither 
numerically or analytically, is much harder. Most stud- 
ies of entanglement in a stationary nonequilibrium set- 
ting have so far considered only small systems of 2 or 3 
spin-1/2 particles [4]. It has been shown that currents 
can enhance entanglement [5j. In the present work we 
present fully analytic results for entanglement in a sta- 
tionary state of a nonequilibrium system of arbitrary size. 
The main questions we study are the dependence of en- 
tanglement in a nonequilibrium stationary state (NESS) 
on the coupling strength to the bath, the driving strength 
measuring how far we are from equilibrium, and the de- 
pendence on system size. The Hamiltonian part of the 
model we consider is a one-dimensional quantum spin 1/2 
chain with XX-type nearest-neighbor coupling, 
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The XX model in its unitary setting is rather simple; 
it can be mapped to a system of noninteracting spinless 
fermions and can be easily diagonalized. A nonequilib- 
rium variant, in which we couple the system to reservoirs 
at chain ends, can also be solved analytically, even in the 
presence of dissipative dephasing. Nonequilibrium dy- 
namics of the XX model with or without dephasing 6 - 
lllj is rather rich, with transport varying from ballistic 
to diffusive. 

For the XX model without dephasing, which shows bal- 
listic transport and is coherent in the bulk, our findings 
are summarized in Figure [l] As one turns on the driv- 
ing fi one moves away from an equilibrium state at infi- 
nite temperature (indicated by an empty circle in Fig.JT]), 
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FIG. 1. (Color online) Expectation values of (cr|) = (z) and 
(<t|ctJ +1 ) = {zz) in the bulk of a NESS for the XX model 
without dephasing. The thick top curve is equilibrium states 
(all at infinite temperature), the shaded region is NESS states 
which can be reached by nonzero driving fi. As one increases 
(j, the NESS state moves downward, eventually reaching a 
region of NESSs having entanglement (dark red region with 
full circle). For details see text. 



eventually reaching NESS in which one has entanglement 
between nearest neighbor spins, independent of system 
size. By a sufficiently strong driving one can induce en- 
tanglement even in an infinite temperature equilibrium 
state. A nonequilibrium setting therefore enlarges the 
region of parameters for which entanglement is present. 
In the presence of dephasing though, situation changes. 
Entanglement is present only if dephasing is sufficiently 
small, or if the system is small enough. In this diffusive 
regime there is no entanglement in the thermodynamic 
limit. 



II. THERMAL ENTANGLEMENT 

Before going to nonequilibrium properties let us first 
remind ourselves of the entanglement in equilibrium 
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state. Results in this section are not new and will serve 
as a reference point against which we will compare NESS 
results. For a more detailed study of entanglement in the 
XY model see, for instance, Ref. [T2]. As a separability 
criterion we will always use the minimal eigenvalue of 
the partially transposed two-spin reduced density matrix 
pJJ+i, where transposition is done with respect to one 
spin. Its negativity is a necessary and sufficient condition 
for entanglement of two spin- 1/2 particles [T3] which is 
the type of entanglement that we study throughout the 
paper. 

We shall consider the grand-canonical state pa at a 
certain inverse temperature j3 = l/T and chemical po- 
tential <j>, 



PG = 



exp (-(3H -4> M) 
Z : 



Z = tr (e 



(-PH-cj>M) 



), (2) 



where M = X^=i a j i s total magnetization. To calcu- 
late the reduced density matrix of two nearest neighbor 
spins in an infinite grand-canonical chain p G we need all 
two spin grand-canonical expectation values (<7"(t" +1 )g, 
where indices a, a' label three Pauli matrices and an 
identity, and (A)q = tr (pqA). Due to the symmetry 
of the XX model the only nonzero terms are (c|)g an d 
(ofaf +l ) G = (°]°\ + i)g - (H) G /(2n) and (a'^ +1 ) G . 
For an infinite chain the energy density and magneti- 
zation can be found in [14) . whereas correlations are in 
Ref. [TS]. In an infinite chain the relation (o-|o-| +1 )g = 
TjCTj+liG holds [15]. In the thermodynamic 
oo we have average magnetization and energy 



limit n - 
density 



i r 

- / tanh (0-2/3 cos &)dfc, (3) 

7T Jo 



cos (k) tanh (<f> — 2(3 cos k) Ak. 



Using these one can easily write a reduced density matrix 
of two nearest-neighbor spins, Pjj+i = \ [1 + (<t|)g(c| ± 

^ + i)+K^i)g(^ + i+<^i)+(^j + i)g^ + i], 
from which one can express the minimal eigenvalue of the 
partially transposed reduced density matrix pJJ+i- The 
result is 



1 



1- 



>G 



(^ z >G-2 



(4) 




FIG. 2. Minimal eigenvalue of pJJ+i Q for a grand-canonical 
state of XX model at <j> = Q. A two-spin reduced density 
matrix Pjj+i is entangled for T < 1.94. 
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FIG. 3. (Color online) Isocurves of constant f3 (solid curves) 
and (j> (dashed curves) in a grand-canonical state pg (pi ■ Ex- 
pectation values of energy density (E) = 2{a*a* + i)G an d 
(z) = (<t|)g |3| in a two-spin nearest neighbor reduced density 
matrix Pj,j+i uniquely determine ft and <j>. Iso-temperature 
curves are at T = 1//3 = 0, 0.5, 1, 2, 4, 8, cx3 (bottom to top), 
while iso-0 curves are at tf> = 0, ±0.25, ±0.5, ±1, ±2 (center to 
right/left). Thick (red) line is expectation values for NESS 
states studied in the present work. 



where w = (o"^cr^ +1 ) G ^ 



Transition from an entan- 



gled state to a separable state happens when A^ n = 0, 
which gives a condition ((<x*<7* +1 )g + V^) 2 = 1 + (a 



z\2 

'jlG- 



Numerically solving this equation for critical f3 at 
we get P c ~ 0.5162, or temperature T c w 1.94. 
This means that below this temperature nearest-neighbor 
spins in a grand-canonical state of the XX model are en- 
tangled, see Fig. [2] Critical j3 varies with <fi only very 
slightly, for instance, at <j) = 4 it is j3 c ps 0.5146. 

Note that the expectation values of magnetization and 
energy density uniquely determine a two-spin nearest- 
neighbor reduced density matrix if the whole system is 



in a grand-canonical state. The relation between these 
two expectation values and temperature T and chemical 
potential 4> can be inferred from Fig. [3] As we shall see, 
all NESS states studied in the present work have energy 
density equal to zero, and therefore lie on the thick (red) 
line in Fig. [3j that is, they have the same expectation 
values of energy and magnetization as equilibrium grand- 
canonical states at infinite temperature, (3 = 0. Note, 
however, that such NESS states are of course not grand- 
canonical as other expectation values, for instance the 
current or (a*ja*j +1 ) G do not have the same expectation 
values as in the grand-canonical state. 
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III. NONEQUILIBRIUM XX MODEL 
A. Setting 

The nonequilibrium dynamics of the XX model will be 
described by the Lindblad master equation [IS] . 



dt 



p = i[p, H] + C h L ath (p) + £ R ath (p) + £ dc ? h (p). (5) 



Each of the three dissipative terms is expressed in terms 
of Lindblad operators Lk as 



C dl %p)=J2(lLkP,L{} + lL k ,pL{] 



(G) 



B. Equilibrium, p = 

The equilibrium situation is obtained for p — (other 
parameters can be arbitrary). The NESS state is, in this 
case, separable and equal to [9], 



Pness 



n 



(12) 



This is exactly the grand-canonical state ^ with j3 — 
and tanh cf> — —p. All equilibrium states, from which 
NESS will be reached for nonzero p, studied in the 
present work therefore are at infinite temperature. These 
are indicated by a thick (red) line in Fig. [3] 



Dephasing is described by n Lindblad operators, each 
acting only on the j-th spin, and being 



- dcph 



(7) 



The bath dissipator £^ ath acts only on the 1st spin while 
£ R a th acts on the last spin. Each involves two Lindblad 
operators, on the left end 



C. Nonequilibrium and no dephasing, 7 = 

We first focus on the case with 7 = because the main 
entanglement features are already present in this simpler 
XX model without dephasing. Because everything is in- 
dependent of the system size n the results presented in 
this subsection are valid for any n. 

Expectation values in the NESS that we need for a 
two-spin reduced density matrix are given in Ref. [5] , 



L\ = X /T L (1 - I + p) *+, L\ = Jr L (l + I - ft) ar 



while on the right end we have 



(8) 



l? = Jr R (i + ^ + p)*+, i* = x /r K (i 



2 



(9) 

where oy = (er* ± i cr J ) / 2 . Relevant parameters are 
the two coupling strengths to the baths Tl,r, driving 
strength p that dictates the magnetization difference be- 
tween chain ends, the average driving p and the dephas- 
ing strength 7. Because driving coefficients have to be 
real we must have Ti^r > as well as (for all four sign 
combinations) 



l±^±p>0. 



(10) 



In order for this to hold p and p must lie inside a rhombus 
with corners (0, ±1) and (±2, 0) in a "p — p" plane. If 
instead of p and p we define the two parameters 



p 



+ p = d, 



2 



p = c, 



(11) 



i.e., p = c + rf, p = (d — c)/2, definition range is a sim- 
ple square, c G [—1,1] and d 6 [—1,1]- The XX model 
with or without dephasing has been studied in Refs. [6j- 
[TTj . Our analytical calculations rely heavily on the exact 
NESS solution found in Refs. [7HS]. 



L ) = t = p 



ai = p 



r L r R 

(i + r L r R )(r L + r R )' 
m (r L -r R ) + r L r R (r L 
2 



(13) 

r R ) 



p 



(i + r L r R )(r L + r R ) 
m (r L -r R )(i-r L r R ) 
2 (i + r L r R )(r L + r R ) 
v (r L -r R )-r L r R (r L + r R ) 
2 (i + r L r R )(r L + r R ) : 



as well as (<TjCr%) — aj(ik~t 2 5j + i^- Ranges that these pa- 
rameters can take are t £ p[0, j], then a\ £ p — ^ [— 1, 1], 
a E p — f [— 1, 1], and a n E p — f [— 1, 1]. 

For equal coupling strengths at both ends, Tl = T R , 
which we will use in all figures, the above expressions 
simplify to t = f , o-i = p — Ti, a = p,, a n = p + Ft. 



1. Two nearest-neighbor spins at the boundary 

The reduced density matrix for two spins at the chain 
end is 



P12 = 



X 














y 


2it 








-2it 


z 

















(14) 



with x = 1 + a + ai + acii — t 2 , y = 1 + a — cii — acii + t 2 , 
z = 1 — a + ai — aoi +i 2 , v = 1 — a — ai + aa\ — 
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t 2 . Calculating pf^ we get only one eigenvalue that can 
possibly be negative and which is 



1 



A min = -(l + aoi-t 



^(a + ai ) 2 +4t2). (15) 



In order to simplify expressions, from now on in this 
subsection we take the same coupling at both ends, 
IY = Tr = r. This has no qualitatively important con- 
sequences. In Fig. [4] we plot in black those regions of 
driving parameters for which p 12 is entangled, that is, 

d" plane we also 



A^j n is negative. In addition to a "c 
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FIG. 5. (Color online) Vertical axis is F (Pl = Tr), horizontal 
axes are c and d. Parameters within the tw o sy mmetrically 
placed shaded regions, existing for T > P m i n ( 18 1, correspond 
to entangled boundary two spins in the NESS. Figure [4] shows 
a cross-section of this plot at T = 2. 



FIG. 4. Parameters for which p\2 is entangled in the NESS 
(black areas). Left: values of c and d (111 for which there 
is entanglement. Right: values of u = Ft = ((<r^) — (o\))/2 
and p with an entangled p\2- Right plot is just a rescaled 
and rotated left figure. Dashed lines in the right plot denote 
a region of allowed parameter values. Tl = Tr = 2. 

show values of u, defined as u = Tt = ((ofj) — (cf))/2 and 
p, which have physical interpretation as the difference of 
magnetization at chain ends and the magnetization in- 
side the chain (at all spins apart from the 1st and the 
last one). On the "c — d" plot the minimal c, i.e., the 
beginning of the right tongue in the left frame of Fig. |4j 
is at 



8r 2 + 4r 4 _ 2r( ' 1 + r 2 ) V8+T2 
i + 2F 2 (i + r 2 ) ' 



(16) 



The functional form of the bottom boundary of this 
tongue, g?(c, r), extending in c from c m i n to c = 1, is 



d = 1 - (1 + w - c(l + 3r 2 + T 4 ) - y/v) 



(17) 



where w = T(l + T 2 )VS + T 2 and v = (1 + T 2 )[(9 



5c 2 )r 4 + (l + c 2 )r b -4(-l+cM;)-2r 2 (-6-2c 2 



■;)]. For 



r > 1 the tongue extends all the way down to d = — 1, 
that is, d(c = 1,T) = —1. For T < 1 it is smaller, 
eventually disappearing for T < T mm 0.511, where 



r s» 1.08873. Regions of T and c and d for which pi2 is 
entangled can be seen in Fig. [5j It is instructive to plot 
how the NESS state differs from the equilibrium grand- 
canonical state. As mentioned, the expectation value of 
the energy is for our NESS states always zero [8] and 
one can always find an equilibrium grand-canonical state 
having the same expectation value of energy and magne- 
tization, see Fig. [3j The expectation value of (<r|cr| +1 ) in 
the NESS on the other hand differs from the one in equi- 
librium. Therefore, we plot in Fig. [6] a region of allowed 
expectation values of a 7 - and OjOj+i m the bulk (away 
from two boundary spins). Remember that in equilib- 
rium at infinite temperature (where (ct^(T^ +1 )g = 0) one 



has the relation (ct|ct| +1 )g 



r z\2 



J/G (top solid curve in 
Fig. [6]) . One can see that for a sufficiently strong driving 
p, keeping p fixed, a NESS is reached in which p 12 is 
entangled (red region). What is even more interesting, 
for sufficiently large (z), that is p, an entangled state 
is reached already for very small values of p. That is, 
for large average magnetization p, one can reach entan- 
glement from an infinite temperature equilibrium state 
already with small driving. One can say that a small 
stationary perturbation of an infinite temperature equi- 
librium state can induce entanglement. This is quite dif- 
ferent from thermal entanglement, which is possible only 
at sufficiently low temperatures. 



3 + (27 - 3V78) 1/3 + (27 + 3V78) 1 / 
6 



(18) 



In other words, P12 can be entangled only if the cou- 
pling strength is in the range r m j n < T < 00. Entangle- 
ment tongues are the largest, i.e., c m i n is the smallest, at 



2. Nearest-neighbor spins in the bulk 

We proceed to entanglement of two nearest-neighbor 
spins in the bulk of the chain. Using again expectation 
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FIG. 6. (Color online) Expectation values of (zz) = (<tJ<tJ +1 ) 
and (z) = {(jj) in the bulk of NESS lie between top solid 
curve and bottom dashed curve. Top solid curve corresponds 
to equilibrium states. As one increases /i away from zero 
(keeping all other parameters fixed) NESS moves vertically 
in the plot. For sufficiently large fi an "entanglement smile" 
region of entangled boundary two spins is reached, denoted 
by the dark (red) region. All is for T = 1.088. 



values ( 14 1 we get the reduced density matrix 




where A± = (1 ± a) 2 - t 2 and B = 1 



(19) 



t 2 . The 



minimal eigenvalue of the partially transposed Pjj+i 
this time 



is 



XH = ±(l + a 2 -t 2 - 2v^T^). (20) 



From now on we again set Tl = Tr = T. In Fig. [7]we can 
see regions of parameters for which Pjj+i is entangled. 
The main difference from the case of pi2 (e.g., Fig. [1]) is 
that the entanglement is present in two pockets in the 
corners of Fig. [7| The lower edge of the right pocket 
(near corner c = d = 1) is 



c(i + 3r 2 + r 4 ) + 2(1 + r 2 )[-rV2 + Va] 



i + r 2 + r 4 



(21) 



where A=l + T{'ST + c 2 T 
from c m ; n to c = 1, with 



r 3 - V8c(l + T 2 )). It runs 

4r(V2 + r(-i + rV2)) 



l + r 2 + r 4 



(22) 



c min is minimal, i.e., the pocket is the largest for V — 
1. When one increases or decreases T the size of this 
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FIG. 7. Regions of entanglement in pj,j+i for Tl = Tr = 1 
in "c — d" (left) and u u — p," plane (right). There are two 
symmetric pockets at large u, that is, at large magnetization 
difference across the chain. 



pocket gets smaller, eventually disappearing when c n 
1. This bounds the range of possible couplings to 



^(1+^-^-1 + 2^2) < T < i(l+V2-^/-l + 2V2). 

(23) 

In the bulk we therefore have entanglement in NESS 
only in this range of T. This can be seen in Fig. [HJ In 




FIG. 8. (Color online) Vertical axis is F (Fl = Tr), horizontal 
are c and d. Two pockets at intermediate couplings 0.531 < 
r < 1.883 (231 are NESSs with entangled two-spin reduced 



density matrix Pt,»+i. 

Fig. [9] we also plot expectation values of magnetization 
and (er|<7| +1 ) in the bulk. One can see, that in contrast 
to entanglement in the boundary two spins, here one has 
entanglement only for sufficiently strong fi. Note how- 
ever that all these states still have zero energy density. 
Compare also the left plot in Fig. [9] with Fig. [T] where 
an enlarged central section of this plot (r = 1) is shown. 

Why can entanglement at the edge of the chain be 
reached already for small /i while large p is needed in the 
bulk? The answer lies in different magnetization values 
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FIG. 9. (Color online) Expectation values in equilibrium, 
p — (top full curve), and out of equilibrium (p ^ 0, region 
between dashed and full curve). Dark (red) region are NESS 
states with bulk nearest-neighbor entanglement. They are 
reached only for sufficiently strong driving p. All is for Tl = 
T R = T, on the left for T = 1, on the right for T = 1.5. An 
enlarged central part of the left frame is shown in Fig. [l] 



at the edge. Because the system is a ballistic transporter 
magnetization exhibits jump at the boundary; it is p, — u 
at the 1st spin and p, at the 2nd, causing entanglement 
for small p. 

Let us finally conclude the section on the XX model 
with a discussion of non nearest-neighbor spins. Calcu- 
lating the reduced density matrix and the minimal eigen- 
value of the partially transposed matrix we see that the 
condition for entanglement is 1 — a 2 < 2y 'a 2 + i 4 , in other 
words, 2\t\ + p 2 > 1. This can never be fulfilled and there 
is no entanglement. 



D. XX with nonzero dephasing 




FIG. 10. Shaded region is a range of coupling strengths F 
(Eq. |25[ ) for which one can have a NESS with nearest-neighbor 
entanglement in the middle of the XX chain with dephasing. 

For the XX model with dephasing NESS expectation 



values depend on the system size. Under nonzero driving 
p the spin current is proportional to ~ similarly as 
long-range correlations. Because dephasing diminishes 
off-diagonal matrix elements in the Pauli basis, it is to be 
expected that the entanglement will be smaller compared 
to a coherent XX model. This is indeed what happens. 

For simplicity we will always use IY = = T and 
look at the central two spins described by the reduced 
density matrix p n /2,n/2+i- Using the exact solution for 
NESS from Refs. 19] one can easily obtain the reduced 
density matrix and the minimal eigenvalue of p^J 2 n/2+v 
Expressions for X^J n as well as boundaries of entangled 
regions are simple but long and we do not give them here. 
The main difference with the XX model is that the size 
of the entanglement pockets depends on the dephasing 
strength 7. There are two symmetrically placed pockets, 
similarly as in the case of no dephasing. Regions of en- 
tangled NESS disappear if 7 is larger than the critical 



7c 



7c 



V2 + 1 



i/r 



(24) 



Also, regions of entangled NESS states are present only 
if the coupling strength T satisfies T_ < T < T + , where 



r+ = 



r (w - j(n - 2) ± y/(w - 7 (n - 2)) 2 - 4) 



(25) 



with w = 1 + \[2. These boundaries, which are a gener- 
alization of Eq.(23), can be seen in Fig. 
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When one increases 7 from zero the last entangled 
NESS to survive at 7 C is the one with c = d = 1, i.e., 
p = 2 and p — 0. At fixed 7 and in the limit n — > 00 
the critical dephasing strength goes to zero and entan- 
glement disappears! With 7 there is therefore a discon- 
tinuous transition from entangled NESSs at 7 = (red 
pocket in Fig. [T]) to no entanglement for 7 > 0. This is 
illustrated in Fig. [TT] One can see how the size of the 
parameter region with entangled NESSs decreases when 
one increases 7 or n. If nearest-neighbor spins Pjj+i are 
in the bulk but not in the middle of chain, j 7^ ra/2, the 



situation is similar as for p. 



n/2,Ti/2+l- 



As one increases 



dephasing 7 or n there is a transition from entangled to 
separable Pj.j+i at certain critical dephasing. The only 
difference is that this critical dephasing weakly depends 
on the location j in the chain where the two spins are. 
Critical j c weakly increases as one moves away from the 
center of the chain. 

For nearest-neighbor spins at the boundary, described 
by the reduced density matrix P12, entanglement is pos- 
sible only if r is larger than some minimal value F m . 
This minimal value depends on 7 and n in a compli- 
cated way, however, it has a very simple thermodynamic 
limit. Namely, linin^oo T m = 1 irrespective of dephasing 
strength 7. For T > 1 dephasing therefore does not de- 
stroy entanglement in pi 2 . The region of entangled states 
still has the same shape as in zero dephasing, Fig. |4j con- 
taining all states with c = 1, only its width c m i n changes. 
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FIG. 11. Left: Dependence of size of entangled region, 
1 — c m i n , on dephasing 7 for three different chain sizes n = 
6, 10, 16. One can see that for 7 > 7 C there are no entangled 
regions (c m i n gets larger than 1). Right: The extent of the 
entangled region for 7 = 0.04 and n = 6 and n = 10 (two full 
curves). c m i n for these two curves is marked with points in 
the left plot. For reference we also plot the entangled region 
without dephasing, 7 = (dashed curve; the same data as in 
Fig. [F|) . The arrow shows the direction in which n increases; 
only one quadrant of the "c — d" plane is shown. Both plots 
are for T = 1 and p n / 2 ,n/2+i- 



It shrinks to zero in the limit n — > 00 as well as for 
7—5-00. In the presence of dephasing there is no entan- 
glement between non- nearest-neighbor spins. 



IV. CONCLUSION 

We have studied bipartite entanglement between two 
spins in a nonequilibrium stationary state of the XX 
chain of any length. In the absence of dephasing we 
show that the two boundary spins can become entan- 
gled for arbitrarily small driving away from an infinite 
temperature equilibrium state. In the bulk, driving must 
surpass a certain value in order to have nearest-neighbor 
entanglement. Coupling with the reservoirs also has to 
be sufficiently strong. Dephasing on the other hand de- 
creases entanglement. Apart from the two boundary 
spins, entanglement is present only if a product of de- 
phasing strength and system size is sufficiently small. 
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